# Analyses for the Web Appendix
library(foreign)
library(readstata13)
library(effects)
library(car)
library(lattice)
library(tidyverse)
library(dplyr)
library(ggthemes)
library(tidyr)
library(plyr)

#### Figures: Predicted probabilities
fill_color <- c("#9cac9d", "#b8d8ba","#cfdbbd","#fcddbc", "#ef959d", "#BF777D","#5F3B3E")

# Figure B3a
## broader equal representation category
df.wide <- read.dta13("coef_wide.dta")
p.HL <- ggplot(data = df.wide) + 
  geom_bar(mapping = aes(x = HL_7_sig_alt1, y = ..prop.., HL_7_sig_alt1 = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#ef959d", "#fcddbc","#cfdbbd" , "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .35),
                     expand = c(0.01, 0)) +
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8),
                     labels=c(expression(beta[H]<~0<~beta[L]), 
                              expression(beta[H]/beta[L]<=1/2), 
                              expression(1/2<~beta[H]/beta[L]<=0.75),
                              expression(0.75<~beta[H]/beta[L]<~1.25), 
                              expression(1.25<=~beta[H]/beta[L]<~2), 
                              expression(2<=beta[H]/beta[L]),
                              expression(beta[L]<~0<~beta[H]),
                              "ambiguous cases")) +
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks = element_blank())
p.HL
ggsave("HL_alt1.pdf",
       width = 250, height = 200, units = "mm")  

# Figure B3b
p.HM <- ggplot(data = df.wide) + 
  geom_bar(mapping = aes(x = HM_7_sig_alt1, y = ..prop.., HM_7_sig_alt1 = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#ef959d", "#fcddbc","#cfdbbd" , "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .35),
                     expand = c(0.01, 0)) +
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8),
                     labels=c(expression(beta[H]<~0<~beta[M]), 
                              expression(beta[H]/beta[M]<=1/2), 
                              expression(1/2<~beta[H]/beta[M]<=0.75),
                              expression(0.75<~beta[H]/beta[M]<~1.25), 
                              expression(1.25<=~beta[H]/beta[M]<~2), 
                              expression(2<=~beta[H]/beta[M]),
                              expression(beta[M]<~0<~beta[H]),
                              "ambiguous cases")) +
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks = element_blank())
p.HM
ggsave("HM_alt1.pdf",
       width = 250, height = 200, units = "mm")  


# Figure B4
### H-L
HL1 <- expression(beta[H]<~0<~beta[L])
HL2 <- expression(beta[H]/beta[L]<=1/2)
HL3 <- expression(1/2<~beta[H]/beta[L]<=0.75)
HL4 <- expression(0.75<~beta[H]/beta[L]<~1.25)
HL5 <- expression(1.25<=~beta[H]/beta[L]<~2)
HL6 <- expression(2<=~beta[H]/beta[L])
HL7 <- expression(beta[L]<~0<~beta[H])

df.all.L <- read.dta13("prpHL_all_alt1.dta")
df.all.L.sorted <- mutate(df.all.L, byvar = reorder(byvar, id))
df.all.L.sorted <- mutate(df.all.L.sorted, re_lab = reorder(re_lab, lab))
df.all.L.sorted$byvar <- revalue(df.all.L.sorted$byvar, c("Issue"="Policy domain"))



p <- ggplot(df.all.L.sorted, aes(x = re_lab, y = re_pr, fill = as.factor(-HL_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HL7, HL6, HL5, HL4, HL3, HL2, HL1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = c(.9, .2),
        legend.justification = c(.9, .2),
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HL_alt1.pdf",
       width = 300, height = 225, units = "mm")  

# Figure B5
### H-M
HM1 <- expression(beta[H]<~0<~beta[M])
HM2 <- expression(beta[H]/beta[M]<=1/2)
HM3 <- expression(1/2<~beta[H]/beta[M]<=0.75)
HM4 <- expression(0.75<~beta[H]/beta[M]<~1.25)
HM5 <- expression(1.25<=~beta[H]/beta[M]<~2)
HM6 <- expression(2<=~beta[H]/beta[M])
HM7 <- expression(beta[M]<~0<~beta[H])

df.all.M <- read.dta13("prpHM_all_alt1.dta")
df.all.M.sorted <- mutate(df.all.M, byvar = reorder(byvar, id))
df.all.M.sorted <- mutate(df.all.M.sorted, re_lab = reorder(re_lab, lab))
df.all.M.sorted$byvar <- revalue(df.all.M.sorted$byvar, c("Issue"="Policy domain"))

p <- ggplot(df.all.M.sorted, aes(x = re_lab, y = re_HMpr, fill = as.factor(-HM_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HM7, HM6, HM5, HM4, HM3, HM2, HM1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = c(.9, .2),
        legend.justification = c(.9, .2),
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HM_alt1.pdf",
       width = 300, height = 225, units = "mm")  




## factor 3 instead of 2
df.wide <- read.dta13("coef_wide.dta")

# Figure B6a
p.HL <- ggplot(data = df.wide) + 
  geom_bar(mapping = aes(x = HL_7_sig_alt2, y = ..prop.., HL_7_sig_alt2 = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#ef959d", "#fcddbc","#cfdbbd" , "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .4),
                     expand = c(0.01, 0)) +
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8),
                     labels=c(expression(beta[H]<~0<~beta[L]), 
                              expression(beta[H]/beta[L]<=1/3), 
                              expression(1/3<~beta[H]/beta[L]<=0.85),
                              expression(0.85<~beta[H]/beta[L]<~1.15), 
                              expression(1.15<=~beta[H]/beta[L]<~3), 
                              expression(3<=beta[H]/beta[L]),
                              expression(beta[L]<~0<~beta[H]),
                              "ambiguous cases")) +
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks = element_blank())
p.HL
ggsave("HL_alt2.pdf",
       width = 250, height = 200, units = "mm")  

# Figure B6b
p.HM <- ggplot(data = df.wide) + 
  geom_bar(mapping = aes(x = HM_7_sig_alt2, y = ..prop.., HM_7_sig_alt2 = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#ef959d", "#fcddbc","#cfdbbd" , "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .4),
                     expand = c(0.01, 0)) +
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8),
                     labels=c(expression(beta[H]<~0<~beta[M]), 
                              expression(beta[H]/beta[M]<=1/3), 
                              expression(1/3<~beta[H]/beta[M]<=0.85),
                              expression(0.85<~beta[H]/beta[M]<~1.15), 
                              expression(1.15<=~beta[H]/beta[M]<~3), 
                              expression(3<=~beta[H]/beta[M]),
                              expression(beta[M]<~0<~beta[H]),
                              "ambiguous cases")) +
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks = element_blank())
p.HM
ggsave("HM_alt2.pdf",
       width = 250, height = 200, units = "mm")  

# Figure B7
### H-L
HL1 <- expression(beta[H]<~0<~beta[L])
HL2 <- expression(beta[H]/beta[L]<=1/3)
HL3 <- expression(1/3<~beta[H]/beta[L]<=0.85)
HL4 <- expression(0.85<~beta[H]/beta[L]<~1.15)
HL5 <- expression(1.15<=~beta[H]/beta[L]<~3)
HL6 <- expression(3<=~beta[H]/beta[L])
HL7 <- expression(beta[L]<~0<~beta[H])

df.all.L <- read.dta13("prpHL_all_alt2.dta")
df.all.L.sorted <- mutate(df.all.L, byvar = reorder(byvar, id))
df.all.L.sorted <- mutate(df.all.L.sorted, re_lab = reorder(re_lab, lab))
df.all.L.sorted$byvar <- revalue(df.all.L.sorted$byvar, c("Issue"="Policy domain"))

p <- ggplot(df.all.L.sorted, aes(x = re_lab, y = re_pr, fill = as.factor(-HL_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HL7, HL6, HL5, HL4, HL3, HL2, HL1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = c(.9, .2),
        legend.justification = c(.9, .2),
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HL_alt2.pdf",
       width = 300, height = 225, units = "mm")  

# Figure B8
### H-M
HM1 <- expression(beta[H]<~0<~beta[M])
HM2 <- expression(beta[H]/beta[M]<=1/3)
HM3 <- expression(1/3<~beta[H]/beta[M]<=0.85)
HM4 <- expression(0.85<~beta[H]/beta[M]<~1.15)
HM5 <- expression(1.15<=~beta[H]/beta[M]<~3)
HM6 <- expression(3<=~beta[H]/beta[M])
HM7 <- expression(beta[M]<~0<~beta[H])

df.all.M <- read.dta13("prpHM_all_alt2.dta")
df.all.M.sorted <- mutate(df.all.M, byvar = reorder(byvar, id))
df.all.M.sorted <- mutate(df.all.M.sorted, re_lab = reorder(re_lab, lab))
df.all.M.sorted$byvar <- revalue(df.all.M.sorted$byvar, c("Issue"="Policy domain"))

p <- ggplot(df.all.M.sorted, aes(x = re_lab, y = re_HMpr, fill = as.factor(-HM_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HM7, HM6, HM5, HM4, HM3, HM2, HM1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = c(.9, .2),
        legend.justification = c(.9, .2),
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HM_alt2.pdf",
       width = 300, height = 225, units = "mm")  



## Broader equal representation category and factor 3 instead of 2
df.wide <- read.dta13("coef_wide.dta")

# Figure B9a
p.HL <- ggplot(data = df.wide) + 
  geom_bar(mapping = aes(x = HL_7_sig_alt3, y = ..prop.., HL_7_sig_alt3 = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#ef959d", "#fcddbc","#cfdbbd" , "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .35),
                     expand = c(0.01, 0)) +
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8),
                     labels=c(expression(beta[H]<~0<~beta[L]), 
                              expression(beta[H]/beta[L]<=1/3), 
                              expression(1/3<~beta[H]/beta[L]<=0.75),
                              expression(0.75<~beta[H]/beta[L]<~1.25), 
                              expression(1.25<=~beta[H]/beta[L]<~3), 
                              expression(3<=beta[H]/beta[L]),
                              expression(beta[L]<~0<~beta[H]),
                              "ambiguous cases")) +
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks = element_blank())
p.HL
ggsave("HL_alt3.pdf",
       width = 250, height = 200, units = "mm")  

# Figure B9b
p.HM <- ggplot(data = df.wide) + 
  geom_bar(mapping = aes(x = HM_7_sig_alt3, y = ..prop.., HM_7_sig_alt3 = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#ef959d", "#fcddbc","#cfdbbd" , "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .35),
                     expand = c(0.01, 0)) +
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8),
                     labels=c(expression(beta[H]<~0<~beta[M]), 
                              expression(beta[H]/beta[M]<=1/3), 
                              expression(1/3<~beta[H]/beta[M]<=0.75),
                              expression(0.75<~beta[H]/beta[M]<~1.25), 
                              expression(1.25<=~beta[H]/beta[M]<~3), 
                              expression(3<=~beta[H]/beta[M]),
                              expression(beta[M]<~0<~beta[H]),
                              "ambiguous cases")) +
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks = element_blank())
p.HM
ggsave("HM_alt3.pdf",
       width = 250, height = 200, units = "mm")  

# Figure B10
### H-L
HL1 <- expression(beta[H]<~0<~beta[L])
HL2 <- expression(beta[H]/beta[L]<=1/3)
HL3 <- expression(1/3<~beta[H]/beta[L]<=0.75)
HL4 <- expression(0.75<~beta[H]/beta[L]<~1.25)
HL5 <- expression(1.25<=~beta[H]/beta[L]<~3)
HL6 <- expression(3<=~beta[H]/beta[L])
HL7 <- expression(beta[L]<~0<~beta[H])

df.all.L <- read.dta13("prpHL_all_alt3.dta")
df.all.L.sorted <- mutate(df.all.L, byvar = reorder(byvar, id))
df.all.L.sorted <- mutate(df.all.L.sorted, re_lab = reorder(re_lab, lab))
df.all.L.sorted$byvar <- revalue(df.all.L.sorted$byvar, c("Issue"="Policy domain"))

p <- ggplot(df.all.L.sorted, aes(x = re_lab, y = re_pr, fill = as.factor(-HL_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HL7, HL6, HL5, HL4, HL3, HL2, HL1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = c(.9, .2),
        legend.justification = c(.9, .2),
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HL_alt3.pdf",
       width = 300, height = 225, units = "mm")  

# Figure B11
### H-M
HM1 <- expression(beta[H]<~0<~beta[M])
HM2 <- expression(beta[H]/beta[M]<=1/3)
HM3 <- expression(1/3<~beta[H]/beta[M]<=0.75)
HM4 <- expression(0.75<~beta[H]/beta[M]<~1.25)
HM5 <- expression(1.25<=~beta[H]/beta[M]<~3)
HM6 <- expression(3<=~beta[H]/beta[M])
HM7 <- expression(beta[M]<~0<~beta[H])

df.all.M <- read.dta13("prpHM_all_alt3.dta")
df.all.M.sorted <- mutate(df.all.M, byvar = reorder(byvar, id))
df.all.M.sorted <- mutate(df.all.M.sorted, re_lab = reorder(re_lab, lab))
df.all.M.sorted$byvar <- revalue(df.all.M.sorted$byvar, c("Issue"="Policy domain"))

p <- ggplot(df.all.M.sorted, aes(x = re_lab, y = re_HMpr, fill = as.factor(-HM_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HM7, HM6, HM5, HM4, HM3, HM2, HM1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = c(.9, .2),
        legend.justification = c(.9, .2),
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HM_alt3.pdf",
       width = 300, height = 225, units = "mm")  

##########################################################################################################################################################
##########################################################################################################################################################
## Five-cat DV
df.wide <- read.dta13("coef_wide.dta")

# Figure B12a
p.HL <- ggplot(data = df.wide) + 
  geom_bar(mapping = aes(x = HL_5_sig, y = ..prop.., HL_5_sig = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#fcddbc", "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .55),
                     expand = c(0.01, 0)) +
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7),
                     labels=c(expression(beta[H]<~0<~beta[L]), 
                              expression(beta[H]/beta[L]<=.85), 
                              expression(0.85<~beta[H]/beta[L]<~1.15), 
                              expression(1.15<=beta[H]/beta[L]),
                              expression(beta[L]<~0<~beta[H]),
                              "", 
                              "ambiguous cases")) +
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks = element_blank())
p.HL
ggsave("HL_5.pdf",
       width = 250, height = 200, units = "mm")  

# Figure B12b
p.HM <- ggplot(data = df.wide) + 
  geom_bar(mapping = aes(x = HM_5_sig, y = ..prop.., HM_5_sig = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#fcddbc", "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .55),
                     expand = c(0.01, 0)) +
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7),
                     labels=c(expression(beta[H]<~0<~beta[M]), 
                              expression(beta[H]/beta[M]<=.85), 
                              expression(0.85<~beta[H]/beta[M]<~1.15), 
                              expression(1.15<=~beta[H]/beta[M]),
                              expression(beta[M]<~0<~beta[H]),
                              "", 
                              "ambiguous cases")) +
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks = element_blank())
p.HM
ggsave("HM_5.pdf",
       width = 250, height = 200, units = "mm")  

## alt color combination
fill_color_5 <- c("#9cac9d", "#b8d8ba","#fcddbc", "#BF777D","#5F3B3E")

# Figure B13
### H-L
HL1 <- expression(beta[H]<~0<~beta[L])
HL2 <- expression(beta[H]/beta[L]<=.85)
HL3 <- expression(0.85<~beta[H]/beta[L]<~1.15)
HL4 <- expression(1.15<=~beta[H]/beta[L])
HL5 <- expression(beta[L]<~0<~beta[H])

df.all.L <- read.dta13("prpHL_all_DV5.dta")
df.all.L.sorted <- mutate(df.all.L, byvar = reorder(byvar, id))
df.all.L.sorted <- mutate(df.all.L.sorted, re_lab = reorder(re_lab, lab))
df.all.L.sorted$byvar <- revalue(df.all.L.sorted$byvar, c("Issue"="Policy domain"))

p <- ggplot(df.all.L.sorted, aes(x = re_lab, y = re_pr, fill = as.factor(-HL_5))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color_5, guide = guide_legend(reverse=FALSE, nrow = 5), name = "", labels = c(HL5, HL4, HL3, HL2, HL1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = c(.9, .2),
        legend.justification = c(.9, .2),
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HL_DV5.pdf",
       width = 300, height = 225, units = "mm")  

# Figure B14
### H-M
HM1 <- expression(beta[H]<~0<~beta[M])
HM2 <- expression(beta[H]/beta[M]<=0.85)
HM3 <- expression(0.85<~beta[H]/beta[M]<~1.15)
HM4 <- expression(1.15<=~beta[H]/beta[M])
HM5 <- expression(beta[M]<~0<~beta[H])

df.all.M <- read.dta13("prpHM_all_DV5.dta")
df.all.M.sorted <- mutate(df.all.M, byvar = reorder(byvar, id))
df.all.M.sorted <- mutate(df.all.M.sorted, re_lab = reorder(re_lab, lab))
df.all.M.sorted$byvar <- revalue(df.all.M.sorted$byvar, c("Issue"="Policy domain"))

p <- ggplot(df.all.M.sorted, aes(x = re_lab, y = re_HMpr, fill = as.factor(-HM_5))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color_5, guide = guide_legend(reverse=FALSE, nrow = 5), name = "", labels = c(HM5, HM4, HM3, HM2, HM1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = c(.9, .2),
        legend.justification = c(.9, .2),
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HM_DV5.pdf",
       width = 300, height = 225, units = "mm")  

######################################################################################################################################################
## Five-cat DV with broader equal representation category
df.wide <- read.dta13("coef_wide.dta")

# Figure B15a
p.HL <- ggplot(data = df.wide) + 
  geom_bar(mapping = aes(x = HL_5_sig_alt1, y = ..prop.., HL_5_sig_alt1 = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#fcddbc", "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .46),
                     expand = c(0.01, 0)) +
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7),
                     labels=c(expression(beta[H]<~0<~beta[L]), 
                              expression(beta[H]/beta[L]<=.85), 
                              expression(0.85<~beta[H]/beta[L]<~1.15), 
                              expression(1.15<=beta[H]/beta[L]),
                              expression(beta[L]<~0<~beta[H]),
                              "", 
                              "ambiguous cases")) +
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks = element_blank())
p.HL
ggsave("HL_5_alt1.pdf",
       width = 250, height = 200, units = "mm")  

# Figure B15b
p.HM <- ggplot(data = df.wide) + 
  geom_bar(mapping = aes(x = HM_5_sig_alt1, y = ..prop.., HM_5_sig_alt1 = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#fcddbc", "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .46),
                     expand = c(0.01, 0)) +
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7),
                     labels=c(expression(beta[H]<~0<~beta[M]), 
                              expression(beta[H]/beta[M]<=.85), 
                              expression(0.85<~beta[H]/beta[M]<~1.15), 
                              expression(1.15<=~beta[H]/beta[M]),
                              expression(beta[M]<~0<~beta[H]),
                              "", 
                              "ambiguous cases")) +
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks = element_blank())
p.HM
ggsave("HM_5_alt1.pdf",
       width = 250, height = 200, units = "mm")  

## alt color combination
fill_color_5 <- c("#9cac9d", "#b8d8ba","#fcddbc", "#BF777D","#5F3B3E")

# Figure B16
### H-L
HL1 <- expression(beta[H]<~0<~beta[L])
HL2 <- expression(beta[H]/beta[L]<=.75)
HL3 <- expression(0.75<~beta[H]/beta[L]<~1.25)
HL4 <- expression(1.25<=~beta[H]/beta[L])
HL5 <- expression(beta[L]<~0<~beta[H])

df.all.L <- read.dta13("prpHL_all_DV5_alt1.dta")
df.all.L.sorted <- mutate(df.all.L, byvar = reorder(byvar, id))
df.all.L.sorted <- mutate(df.all.L.sorted, re_lab = reorder(re_lab, lab))
df.all.L.sorted$byvar <- revalue(df.all.L.sorted$byvar, c("Issue"="Policy domain"))

p <- ggplot(df.all.L.sorted, aes(x = re_lab, y = re_pr, fill = as.factor(-HL_5))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color_5, guide = guide_legend(reverse=FALSE, nrow = 5), name = "", labels = c(HL5, HL4, HL3, HL2, HL1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = c(.9, .2),
        legend.justification = c(.9, .2),
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HL_DV5_alt1.pdf",
       width = 300, height = 225, units = "mm")  


# Figure B17
### H-M
HM1 <- expression(beta[H]<~0<~beta[M])
HM2 <- expression(beta[H]/beta[M]<=0.75)
HM3 <- expression(0.75<~beta[H]/beta[M]<~1.25)
HM4 <- expression(1.25<=~beta[H]/beta[M])
HM5 <- expression(beta[M]<~0<~beta[H])

df.all.M <- read.dta13("prpHM_all_DV5_alt1.dta")
df.all.M.sorted <- mutate(df.all.M, byvar = reorder(byvar, id))
df.all.M.sorted <- mutate(df.all.M.sorted, re_lab = reorder(re_lab, lab))
df.all.M.sorted$byvar <- revalue(df.all.M.sorted$byvar, c("Issue"="Policy domain"))


p <- ggplot(df.all.M.sorted, aes(x = re_lab, y = re_HMpr, fill = as.factor(-HM_5))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color_5, guide = guide_legend(reverse=FALSE, nrow = 5), name = "", labels = c(HM5, HM4, HM3, HM2, HM1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = c(.9, .2),
        legend.justification = c(.9, .2),
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HM_DV5_alt1.pdf",
       width = 300, height = 225, units = "mm")  


####################################################################################
## Fixed Effects models

# Figure C1
### H-L
HL1 <- expression(beta[H]<~0<~beta[L])
HL2 <- expression(beta[H]/beta[L]<=1/2)
HL3 <- expression(1/2<~beta[H]/beta[L]<=0.85)
HL4 <- expression(0.85<~beta[H]/beta[L]<~1.15)
HL5 <- expression(1.15<=~beta[H]/beta[L]<~2)
HL6 <- expression(2<=~beta[H]/beta[L])
HL7 <- expression(beta[L]<~0<~beta[H])

df.all.L <- read.dta13("prpHL_all_FE.dta")
df.all.L.sorted <- mutate(df.all.L, byvar = reorder(byvar, id))
df.all.L.sorted <- mutate(df.all.L.sorted, re_lab = reorder(re_lab, lab))
df.all.L.sorted$byvar <- revalue(df.all.L.sorted$byvar, c("Issue"="Policy domain"))

p <- ggplot(df.all.L.sorted, aes(x = re_lab, y = re_pr, fill = as.factor(-HL_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HL7, HL6, HL5, HL4, HL3, HL2, HL1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = "right",
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HL_FE.pdf",
       width = 300, height = 225, units = "mm")  

# Figure C2
### H-M
HM1 <- expression(beta[H]<~0<~beta[M])
HM2 <- expression(beta[H]/beta[M]<=1/2)
HM3 <- expression(1/2<~beta[H]/beta[M]<=0.85)
HM4 <- expression(0.85<~beta[H]/beta[M]<~1.15)
HM5 <- expression(1.15<=~beta[H]/beta[M]<~2)
HM6 <- expression(2<=~beta[H]/beta[M])
HM7 <- expression(beta[M]<~0<~beta[H])

df.all.M <- read.dta13("prpHM_all_FE.dta")
df.all.M.sorted <- mutate(df.all.M, byvar = reorder(byvar, id))
df.all.M.sorted <- mutate(df.all.M.sorted, re_lab = reorder(re_lab, lab))
df.all.M.sorted$byvar <- revalue(df.all.M.sorted$byvar, c("Issue"="Policy domain"))

p <- ggplot(df.all.M.sorted, aes(x = re_lab, y = re_HMpr, fill = as.factor(-HM_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HM7, HM6, HM5, HM4, HM3, HM2, HM1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = "right", 
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HM_FE.pdf",
       width = 300, height = 225, units = "mm")  

##############################################################################################################################
##############################################################################################################################
#### US studies only

# Figure E1
df.wide <- read.dta13("coef_wide.dta")
df.wide2 <- dplyr::filter(df.wide, t_L<10 & t_M<10 & t_H<10  & us==1)

p.dens <- ggplot(data=df.wide2) +
  geom_histogram(aes(x =t_L), bins=50, linetype=0, colour="red", fill="red", alpha=.3) +
  geom_histogram(aes(x =t_M), bins=50, linetype=0, colour="darkgreen", fill="darkgreen", alpha=.3) +
  geom_histogram(aes(x =t_H), bins=50, linetype=0, colour="blue", fill="blue", alpha=.3) +
  scale_x_continuous(breaks = c(-5,-1.96,0,1.96,5,10,15), labels = c("-5","-1.96","0","1.96","5","10","15")) +
  labs(x = expression(paste(beta/se)), y="")+
  theme_bw() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=18))
p.dens + 
  annotate(geom="text", x=-2, y=15.75, label="Low-income test statisics",
           size=7.5,
           color="red", 
           alpha=.9) +
  geom_curve(aes(x=-2.65,y=14.75,xend=-1.5,yend=13),
             curvature = -0.1, linetype="solid", size=.1, color="red",
             arrow = arrow(type="open",length = unit(0.25,"cm"))) +
  annotate(geom="text", x=1.4, y=16.5, label="Middle-income test statisics",
           size=7.5,
           color="darkgreen", 
           alpha=.9) +
  geom_curve(aes(x=1.3,y=15.85,xend=.85,yend=13.25),
             curvature = 0.1, linetype="solid", size=.1, color="darkgreen",
             arrow = arrow(type="open",length = unit(0.25,"cm"))) +
  annotate(geom="text", x=6.25, y=15.75, label="High-income test statisics",
           size=7.5,
           color="blue", 
           alpha=.9) +
  geom_curve(aes(x=3.9,y=14.75,xend=2.9,yend=13),
             curvature = 0.1, linetype="solid", size=.1, color="blue",
             arrow = arrow(type="open",length = unit(0.25,"cm"))) 
ggsave("t_statistics_US.pdf", 
       width = 250, height = 200, units = "mm")


## 
df.wide <- read.dta13("coef_wide.dta")
df.US <- dplyr::filter(df.wide, us==1)

# Figure E2a
p.HL <- ggplot(data = df.US) + 
  geom_bar(mapping = aes(x = HL_7_sig, y = ..prop.., HL_7_sig = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#ef959d", "#fcddbc","#cfdbbd" , "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .35),
                     expand = c(0.01, 0)) +
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8),
                     labels=c(expression(beta[H]<~0<~beta[L]), 
                              expression(beta[H]/beta[L]<=1/2), 
                              expression(1/2<~beta[H]/beta[L]<=0.85),
                              expression(0.85<~beta[H]/beta[L]<~1.15), 
                              expression(1.15<=~beta[H]/beta[L]<~2), 
                              expression(2<=beta[H]/beta[L]),
                              expression(beta[L]<~0<~beta[H]),
                              "ambiguous cases")) +
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks = element_blank())
p.HL
ggsave("HL_US.pdf",
       width = 250, height = 200, units = "mm")  

# Figure E2b
p.HM <- ggplot(data = df.US) + 
  geom_bar(mapping = aes(x = HM_7_sig, y = ..prop.., HM_7_sig = 1), 
           stat = "count", width = .9, color="black", 
           fill=c("#5F3B3E", "#BF777D", "#ef959d", "#fcddbc","#cfdbbd" , "#b8d8ba", "#9cac9d", "#a59cac"),
           size=1, 
           alpha=1)+
  ylab("Proportion") +
  xlab("") +
  scale_y_continuous(limits = c(0, .35),
                     expand = c(0.01, 0)) +
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8),
                     labels=c(expression(beta[H]<~0<~beta[M]), 
                              expression(beta[H]/beta[M]<=1/2), 
                              expression(1/2<~beta[H]/beta[M]<=0.85),
                              expression(0.85<~beta[H]/beta[M]<~1.15), 
                              expression(1.15<=~beta[H]/beta[M]<~2), 
                              expression(2<=~beta[H]/beta[M]),
                              expression(beta[M]<~0<~beta[H]),
                              "ambiguous cases")) +
  theme_bw() +
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks = element_blank())
p.HM
ggsave("HM_US.pdf",
       width = 250, height = 200, units = "mm")  


# Figure E3
### H-L
HL1 <- expression(beta[H]<~0<~beta[L])
HL2 <- expression(beta[H]/beta[L]<=1/2)
HL3 <- expression(1/2<~beta[H]/beta[L]<=0.85)
HL4 <- expression(0.85<~beta[H]/beta[L]<~1.15)
HL5 <- expression(1.15<=~beta[H]/beta[L]<~2)
HL6 <- expression(2<=~beta[H]/beta[L])
HL7 <- expression(beta[L]<~0<~beta[H])

df.all.L <- read.dta13("prpHL_all_US.dta")
df.all.L.sorted <- mutate(df.all.L, byvar = reorder(byvar, id))
df.all.L.sorted <- mutate(df.all.L.sorted, re_lab = reorder(re_lab, lab))
df.all.L.sorted$byvar <- revalue(df.all.L.sorted$byvar, c("Issue"="Policy domain"))


p <- ggplot(df.all.L.sorted, aes(x = re_lab, y = re_pr, fill = as.factor(-HL_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HL7, HL6, HL5, HL4, HL3, HL2, HL1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = "right",
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HL_US.pdf",
       width = 300, height = 225, units = "mm")  


# Figure E4
### H-M
HM1 <- expression(beta[H]<~0<~beta[M])
HM2 <- expression(beta[H]/beta[M]<=1/2)
HM3 <- expression(1/2<~beta[H]/beta[M]<=0.85)
HM4 <- expression(0.85<~beta[H]/beta[M]<~1.15)
HM5 <- expression(1.15<=~beta[H]/beta[M]<~2)
HM6 <- expression(2<=~beta[H]/beta[M])
HM7 <- expression(beta[M]<~0<~beta[H])

df.all.M <- read.dta13("prpHM_all_US.dta")
df.all.M.sorted <- mutate(df.all.M, byvar = reorder(byvar, id))
df.all.M.sorted <- mutate(df.all.M.sorted, re_lab = reorder(re_lab, lab))
df.all.M.sorted$byvar <- revalue(df.all.M.sorted$byvar, c("Issue"="Policy domain"))

p <- ggplot(df.all.M.sorted, aes(x = re_lab, y = re_HMpr, fill = as.factor(-HM_7))) + 
  geom_bar(position="fill", stat = "identity")+ 
  facet_wrap(.~byvar, scales="free", nrow=2)+
  scale_fill_manual(values= fill_color, guide = guide_legend(reverse=FALSE, nrow = 7), name = "", labels = c(HM7, HM6, HM5, HM4, HM3, HM2, HM1)) +
  labs(title="", x = "", y="Predicted probability", fill="")+
  ylab("Predicted probability") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        aspect.ratio = 1, 
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = "right",
        legend.text = element_text(size=13),
        strip.text.x = element_text(size = 14, face="bold"))
p 

ggsave("all_HM_US.pdf",
       width = 300, height = 225, units = "mm")  
